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Abstract 

A non-isothermal phase field model that captures both displacive and diffusive phase transfor- 
mations in a unified framework is presented. The model is developed in a formal thermodynamic 
setting, which provides guidance on admissible constitutive relationships and on the coupling of the 
numerous physical processes that are active. Phase changes are driven by temperature-dependent 
free-energy functions that become non-convex below a transition temperature. Higher-order spatial 
gradients are present in the model to account for phase boundary energy, and these terms neces- 
sitate the introduction of non-standard terms in the energy balance equation in order to satisfy 
the classical entropy inequality point-wise. To solve the resulting balance equations, a Galerkin 
finite element scheme is elaborated. To deal rigorously with the presence of high-order spatial 
derivatives associated with surface energies at phase boundaries in both the momentum and mass 
balance equations, some novel numerical approaches are used. Numerical examples are presented 
that consider boundary cooling of a domain at different rates, and these results demonstrate that 
the model can qualitatively reproduce the evolution of microstructural features that are observed 
in some alloys, especially steels. The proposed model opens a number of interesting possibilities 
for simulating and controlling microstructure pattern development under combinations of thermal 
and mechanical loading. 

Keywords: Displacive transformations, diffusive transformations, martensite, pearlite, 
thermodynamics, phase field models, finite element methods. 



1. Introduction 

Displacive and diffusive phase transformations in solids, driven by temperature changes or me- 
chanical loading, occur in many industrial and natural processes. In particular, careful temperature 
control is used extensively to architect microstructural features of alloys and thereby tailor their 
mechanical properties. We present in this work a non-isothermal phase field model for simulating 
the development of microstructure due to both displacive and diffusive processes under thermal 
and mechanical loading. Continuum-level order parameters are defined that indicate the presence 
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of various phases, and coupled evolution equations for these parameters are developed. The model 
is framed in a thermodynamic setting which provides guidance on the coupling of various processes. 



Displacive phase transitions have been extensively studied over a long period. Wechsler et al. 



( 1953 ) and Bowles and McKenzie ( 1954 ) proposed theories of cubic-to-tetragonal martensitic phase 



transitions in which the key feature is the description of the 'Bain distortion' (Bhadeshia 1987 



Wayman 


1990). The formation of martensite twins has also been studied 


minimisers ( 


Ball and James 


1987 


Kohn 


1991, 


Bhattacharya 


1991 


). Falk ( 



proposed a model 

based on a single order parameter, with a free-energy that is a non-convex sextic polynomial in the 
order parameter, although the model did not include evolution equations for the order parameter 
or surface (phase boundary) energy contributions to the free-energy. It is not uncommon that such 
singular models are studied, see for example Ball (2004). However, models that do not account 
for surface energy cannot predict the detail of the microstructure, which is determined by the 
relative balance between bulk and interface energy, and such models will generally lead to an ill- 
posed problem when inserted into a differential balance equation. The above references do not 



consider evolution equations, in which case kinetic aspects of a transition are ignored. Barsch 



and Krumhansl (1984) and Jacobs (1985) proposed models featuring both bulk and interfacial 



contributions to the free-energy, together with balance of linear momentum, thereby addressing 
the evolution of transitions. Our treatment of the displacive transformations is in the same spirit 
as Barsch and Krumhansl (1984) and Jacobs (1985). There exist other phase field approaches to 



martensitic transformations that include bulk and surface energy, but that do not invoke balance 
of momentum (Wang and Khachaturyan 1997 Shenoy et al. 1999 Onuki 1999). These models 
involve the solution of a diffusion-type problem. 

For diffusive transformations, the phase field approach has been used extensively to model 
microstructure evolution. The best known model for conserved order parameters is the Cahn- 



Hilliard equation (Cahn and Hilliard 1958), in which surface energy is introduced via spatial 



gradients of the order parameter in the free-energy. Of relevance to the problems that we will 
consider, the Cahn-Hilliard equations was extended to include non-isothermal effects in Alt and 



Pawlow (1992). 



In the heat treatment of steel alloys, both displacive and diffusive transformation are im- 
portant. Displacive transformations (martensite formation) occur rapidly compared to diffusive 
transformations (pearlite formation), with the speed of latter being limited by the rate at which 
carbon atoms can diffuse. There have been many studies into models for displacive or diffusive 
phase transitions, but fewer into models that can capture both displacive and diffusive transforma- 
tions, and interactions between the resulting phases. Wang et al. (1993) presented an isothermal 
model Ginzburg-Landau type model for simulating both displacive and diffusive transformations 
via coupled Cahn-Allen and Cahn-Hilliard type equations. Consistent with the Ginzburg-Landau 



framework, the model of Wang et al. ( 1993 ) does not invoke balance of momentum. Bouville and 



Ahluwalia (2006) presented an isothermal model for coupled displacive/diffusive processes (see also 



Bouville and Ahluwalia (2007)), in which the treatment of displacive processes resembles that of 



Barsch and Krumhansl (1984), and diffusive processes are modelled by a Cahn-Hilliard type equa- 



tion. Interactions between processes was accounted for via coupling terms in the free-energy, and 
the impact of various coupling terms was studied numerically. 

While the heat treatment of steel is perhaps the most technologically relevant, diffusive and 
displacive transformations driven by mechanical and / or temperature loading are relevant for a va- 
riety of materials. For example, experimental observations of martensitic twinning in perovskites 
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are presented by Harrison et al. (2004). Our intention in this work is to present a generic frame- 
work and to demonstrate that the development of key microstructural features can be modelled. 
Because of the degree to which heat treatment is used for steels, we will at times use steel-specific 
terminology and refer to a high temperature stable phase as austenite and a diffusive phase as 
'pearlite', despite the generic nature of the formulation. 

A focus of this work is the formulation of a non-isothermal model in a thermodynamic setting. 
Free-energy expressions that involve non-standard higher-order gradient terms to model surface 
energies are postulated, and a model that satisfies the first and second laws of thermodynamics is 
formulated. Dealing with the higher-order gradients in a classical thermodynamic setting is not 
trivial, so in formulating our model we adopt an approach with parallels to that advocated by 



Gurtin (1996), and we show that the concept of 'work done' by higher-order (nonlocal) terms on 



the boundary of a domain is necessary to formulate a model that satisfies the entropy inequality 
point-wise. Many of the aforementioned works on phase field models invoke thermodynamic argu- 
ments, but usually premised on the a priori postulation of a chemical potential and proportionality 
between a mass flux and the chemical potential. Following ad-hoc arguments Maraldi et al. (2010) 
include the heat equation in their phase field model, following the usual methodology of assuming 
a structure for the stress tensors and the chemical potential. A distinguishing feature of our for- 
mulation is that it follows from an energy balance that is posed in terms of the work done on the 
boundary of a domain, in the spirit of Gurtin ( 1996 ) . The equations that result from our derivation 
are not trivial to solve. In particular, the incorporation of surface energies via higher-order spatial 
gradients leads to coupled fourth-order hyperbolic and parabolic equations. We therefore pay at- 
tention to the development of a Galerkin finite element formulation of the equations, which we then 
use to compute a number of example problems. The presented formulation is novel in the coupling 
of displacive, diffusive and thermal processes. Moreover it is distinguished from other works by the 
formal thermodynamic framework in which the model is developed, and the sophistication of the 
numerical method formulated for solving the resulting equations. 

The rest of this work is structured as follows. The scale of the problem that we consider 
is discussed and order parameters that are relevant at the considered scale are defined. This is 
followed by the formulation of thermodynamic balance laws and formalisation of the restrictions 
imposed by the entropy inequality. Balance laws and constitutive models are then synthesised to 
yield governing equations, which is followed by the development of a suitable Galerkin finite element 
method. A number of simulations that qualitatively illustrate features of non-isothermal coupled 
phase transformations are presented. The computer code used to produce all numerical examples 
is freely available and distributed under a GNU public license. It can be found in the supporting 
material (Wells and Maraldi . 2011 ). Following the numerical examples, some conclusions are drawn. 



2. Order parameters and scale of the problem 

We define now order parameters whose values will be used to identify phases. The choice of 
order parameters depends on the scale of observation. We consider the level of a single grain, at 
which scale the lattice orientation can be inferred, but at which the continuum hypothesis remains 
valid. The domain of interest will be denoted by C M. d , where d is the spatial dimension. 

Diffusive-type transformations will be characterised by a conserved scalar order parameter c. 
The parameter c typically represents the relative concentration of an alloying element, such as 
carbon, and can act as an indicator of phase. In steels, for example, if c represents the deviation 
in the carbon concentration away the equilibrium concentration in austenite, then the presence of 
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pearlite is indicated by regions in which c alternates spatially about zero. The carbon-rich cementite 
phase can be identified by positive c, and the ferrite phase can be identified by negative c. 

To characterise displacive transformations, deformation-related order parameters are consid- 
ered. In particular, we will use scalar order parameters that are functions of the linearised 
strain tensor e = ( Vtt + (Vm) T ) /2, where u is the displacement vector. The formation of phases 



driven by displacive effects at the scale of observation that we have chosen is dependent on the 
lattice orientation. We restrict ourselves to two spatial dimensions (d = 2) and an initially square 
lattice aligned with the Cartesian x\ and x 2 axes, in which case we consider a volumetric order 
parameter e%, 

ei=en+e 2 2, (1) 

an order parameter e 2 , 

&2 = £11 - &22, (2) 

and the shear strain e^, 

e3 = £12- (3) 
The order parameter e 2 acts as an indicator of a martensitic phase, with a significant variation 



away from zero indicating a square-to- rectangle type transition (Jacobs, 1985). Martensite twins 
will be evident when e 2 alternates spatially about zero. 

3. Thermodynamic balance laws and restrictions 

To formulate thermodynamic restrictions on constitutive equations, it is sufficient at this stage 
to postulate the existence of a suitably smooth Helmholtz free-energy density functional / of the 
form 

/ = /(e,Ve,c,Vc,T), (4) 

where T is the temperature. A dependency of the free-energy on the gradient terms Ve and Vc is 
included with a view to the inclusion of surface energies in the context of a diffuse interface model. 
A precise form of the free-energy as a function of the order parameters introduced in Section [2] will 
be presented in Section [5} 

Special care is required in constructing the energy balance in order to properly account for 
mass diffusion and the presence of higher-order spatial gradients in the free-energy. We adopt an 



approach to the energy balance that shares features with the approach of Gurtin (1996) for the 



Cahn-Hilliard equation, in which non-standard force-like terms associated with higher order terms 



are identified and accounted for in the energy balance equation. Unlike in the work of Gurtin 



(1996), we do not specify a priori any new balance laws for these extra terms, but we will show 
that balance laws for these terms are implied by insisting upon satisfaction of the classical entropy 
inequality. 

We will consider linearised kinematics throughout. The open domain R C f2 is used to denote an 
arbitrary sub-region of fi, hence integral balance laws posed on R can be localised. The outward 
unit normal vector to R (on dR) is denoted by n. Since c is assumed to be a conserved order 
parameter it must satisfy 



4 



where j c is the mass flux of c. The classical linear momentum balance equation reads: 



— / pv dx = f V • cr dx + f b dx, 
dt Jr Jr Jr 



(6) 



where p is the mass density, assumed to be constant, v is the velocity, cr = cr T is the stress and b 
is a body force. Localising equation ([6]), 



dv 



V • cr + b. 



(7) 



Multiplying the balance of linear momentum equation Q by u, and then integrating over R and 
applying integration by parts leads to the mechanical energy balance: 



ld_ 
2dt 



/ p||u|| 2 dx = — / a:Vvdx+ I t ■ v ds + / b-vdx, 
Jr Jr JdR Jr 



(8) 



where t = an is the traction. 

We now consider an energy balance equation of the form 



d 
dt 



• dx 



/ - p\\v\\ 2 + u dx = — / q-nds + / t-vds + / bvi 
Jr 2 JdR JdR Jr 

+ / En: ids + / $,-ncds— pj c -nds, (9) 
JdR JdR JdR 

where u is the specific internal energy density (per unit volume), q is the heat flux, the third-order 
tensor 5] and the vector £ are stress-like terms, and pj c is related to the energy transported into 
the domain by the mass flux of c. The terms E and £ are not standard, and their presence is a 
consequence of the nonlocality implied by the dependency of the free-energy on Vs and Vc. The 
'fluxes' Sn : e and £ • nc, evaluated on dR, can be interpreted as the power expended across the 



boundary of R by a material particle just outside R. This concept is used by Gurtin (1996) for the 



Cahn-Hilliard equation, and less explicitly by Polizzotto (2003) for gradient elasticity problems 



Polizzotto (2003) introduces the concept of a 'nonlocal' residual in the energy balance equation to 



account the non-standard terms that arise in gradient elasticity. The precise role of SI and £ will 
become more evident when considering admissible constitutive equations. Inserting the mechanical 
energy balance ^ into the energy balance ^ and applying the divergence the theorem leads to a 
local internal energy balance equation: 

u = -V • q + cr : Vv + (V • S) : Vv + S : Ve + (V • £)c + £ • Vc - pV ■ j c - Vp ■ j c . (10) 

We will insist upon the satisfaction of a local entropy inequality of the form 

(11) 

where s is the entropy density per unit volume. Using the definition of the Helmholtz free-energy 



s>-V .(J 



(/ = u — Ts) in the entropy inequality ( 11 ) leads to 



u-f-Ts + TV- (J) >(). 



(12) 
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Inserting the internal energy balance (10) into the above inequality yields 



-V-g + tr:V« + V- £:e + £:Ve + V-£c + £-Vc 



+ f id-V f i-j c -f-fs + TV- (|) >0, (13) 



which can be manipulated into the form 



cr : Vv + V • £ : k + £ : Ve + V • £c + £ • Vc + \ic - • j c - / - Ts - VT • (|;) > 0. (14) 
We will insist upon satisfaction of this inequality in our model. 

4. Admissible constitutive equations 



By insisting upon satisfaction of the entropy inequality (14), the precise form of some consti- 
tutive models will follow directly from the definition of the free-energy, while for the others it will 
simply imply a restriction. To start, for a Helmholtz free-energy that has the functional form of 
equation Q, its time derivative reads: 



Inserting the above expansion of / into the entropy inequality (14), 



°° + VS -f^- i + °'- i+ (*-Me)- Vi -{ S+ %) f 

+ ( v .« + (1 -|) 6+ ( f -^ ; ).Vc-VM- J -e-Vr.(|)>0, (16) 

where is has been assumed that the stress tensor can can be decomposed additively into inviscid 
(<r e ) and viscous (<r„) partsQ 

cr = a,, + rr r . (17) 
From the arbitrariness of V, c and T, and the insistence upon satisfaction of equation (1 161) , we 



can infer various admissible constitutive relationships. Equation (16) implies for the 'higher-order' 
stress I] that 

s = Wi> (18) 

for the inviscid component of the stress that 

^-fe-^ 5 -^"^^' (19) 
and that a constitutive model for the viscous stress must satisfy 

<t v : Vv > 0. (20) 



J This permits the incorporation of a Kelvin- Voigt type model. Other models that involve springs and dashpots 
in series can be formulated via the introduction of strain-like internal variables. 
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Satisfaction of equation (16) also requires that the 'chemical higher-order stress' be given by 



*-m- (21) 



that the 'chemical potential' is given by 



and that a constitutive model for the mass flux vector must satisfy 

V/i • jc < 0- (23) 

The entropy density is given by 

M 

For the heat flux, the entropy inequality requires that 

VT-(|)<0. (25) 
Without the inclusion of the non-standard terms S and £ in the energy balance, the classical 



entropy inequality could not be satisfied point-wise, as shown in Gurtin (1965) for elasticity. This 



is evident in equation (16), which in the absence of S and £ could not be guaranteed to hold since 



the terms df/dVe and df/dVc would be without a 'partner' stress term. 

5. Specific form of the Helmholtz free-energy and constitutive equations 

Before presenting the boundary value problems that define the complete model, it is useful 
to specify more precisely the functional form of a Helmholtz free-energy density in terms of the 
order parameters that is suitable for modelling phase transformations. It also useful to define the 
adopted constitutive models that do not follow as a direct consequence of the chosen free-energy 
expression. While the governing equations will be presented in a format that is largely independent 
of the details of the free-energy function, some poignant features of the governing equations only 
become apparent after the introduction of particular constitutive equations, especially those related 
to the surface energy. Where convenient, we will express constitutive models in terms of derivatives 
of the free-energy. This is because the presented numerical simulations employ novel techniques 
for the automated generation of computer code from a domain-specific language that can compute 
the necessary derivatives automatically. The computer code therefore requires expressions for the 
free-energy only. 

5.1. Helmholtz free- energy 

We consider a Helmholtz free-energy density that can be additively composed according to 

/ = /diff (T, c, Vc) + /disp (T, ei, e 2 , e 3 , Ve 2 , c) + f cp \ (c, e 2 ) + /therm (T) , (26) 

where /diff is the free-energy associated with diffusive processes, /disp is the free-energy associated 
with displacive transformations and mechanical deformation, / cp j is the free-energy associated with 
the interaction of phases and /therm is a part of the free-energy which is dependent on the temper- 
ature only. The strain-related order parameters ei and the mass concentration order parameter c 
were defined in Section [2| Th e functional forms of /diff, /disp and / cp i that we adopt come from 



Bouville and Ahluwalia (2006), with some minor modifications. 



xlO 
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Figure 1: Diffusive free-energy density as a function of c at various temperatures T (Aa = 7.31 x 10 3 , A2 
6.62 x 1(T 3 , T P = 1). 

5.1.1. Diffusive part of the free-energy 

We postulate a diffusive free-energy function of the form 



Ai 4 A 2 T-T P 2 A c 
/ diff = T c + T ^^c + y | 



Veil 



(27) 



where A4, A2 and A c are positive constants, and Tp is the non-dimensional temperature above 
which /jiff is convex in c. In the context of steel, Tp is the temperature above which austenite is 
the stable phase. For T < Tp, f^m becomes a double-well function. This can be seen in Figure [TJ 
in which the diffusive free-energy as a function of c is plotted for various temperatures. The further 
the temperature drops below Tp, the deeper the energy wells. The existence of two wells is what 
can lead to regions with layers alternating values of c, which is typical of a pearlitic structure. The 
presence of the term Vc in the free-energy accounts for the energy associated with the formation 
a phase boundaries. In the context of pearlite, it reflects the surface energy associated with the 
boundaries between cementite and ferrite. Mathematically, it provides a regularising effect when 
the free-energy is non-convex with respect to c. Note that below Tp there is no stable local 
minimum at c = for the chosen free-energy function. 

5.1.2. Displacive part of the free-energy 

The displacive part of the free-energy is postulated as follows: 



disp 



Bq2 6 

IT 2 



-B42 4 B22 T 

e% H 

4 2 2 T M 



' M e\ + [ei - (a (T - T ref ) + x\ c c + x^ej)] 



+ 



B. 



3„2 



I Ve 3 | 



(28) 



where Bq2, B42, B22, B±, B% and A e are positive constants, Tm is the non-dimensional temperature 
below which /di sp is non-convex in e2, a?i c and X12 are constant coupling parameters that induce 
volumetric changes as a consequence of diffusive and displacive phase changes, respectively, and 
a > is the thermoelastic coefficient which determines volumetric changes as a consequence of 
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temperature deviations away from a reference temperature T re f. The volumetric term in the dis- 
placive free-energy is chosen such that the model will coincide with classical thermoelasticity. The 
displacive contribution to the free-energy as a function of e 2 at various temperatures is illustrated 
in Figure [2] Physically, Tm is the temperature below which martensite can form. The further 
the temperature drops below Tm, the deeper the double- wells in /d isp as a function of e 2 - This 
is what leads to formation of martensite, with spatially alternating values of e 2 indicating twin- 
ning. Analogous to the diffusive contribution to the free-energy, the gradient term Ve 2 accounts 
for the energy associated with twin boundaries, and provides a regularising effect at temperatures 
below Tm- Similar to the diffusive free-energy contribution, below Tm there is no stable local 
minimum at e2 = for the chosen free-energy function. 

5.1.3. Thermal part of the free- energy 

The thermal part of the free-energy is postulated to be of the classical form 

/therm = —Sq{T — T re f) — — " (T — T re f ) 2 , (29) 

where sq is the reference entropy density and c v is the heat capacity, both of which are assumed 
to be positive and constant, and T re f is the aforementioned constant reference temperature. 

5.1.4- Phase coupling contribution 

Displacive and diffusive phases do not generally co-exists at a given point. To model this, we 
include in the free-energy a contribution of the form 

/cpi = x 2c c 2 e 2 2 , (30) 

where x 2c > is a constant. The task of this term is to penalise energetically concurrent variations 
away from zero of e 2 and c. That is, it penalises the co-existence of pearlitic (c 7^ 0) and martensitic 
phases [e 2 / 0). 
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5.2. Constitutive models as a consequence of the free- energy 

The thermodynamic restrictions in Section [4j together with the Helmholtz free-energy defined 
in this section, provide constitutive models for the stress, the chemical potential and the entropy. 
We provide now some expansions for these constitutive models in terms of the order parameters 
for the particular free-energy that we consider. 

Introducing the notation <r = df/ds, for the 'local' part of the inviscid stress tensor (see 
equation (19)), a in terms of derivatives of / with respect to the order parameters reads: 

df df df 

<x = -w— (ei <g> ei + e 2 <g> e 2 ) + 77— (ei ® e\ — e 2 <g> e 2 ) + 77— e 3 [e.\ ® e 2 + e 2 ® ei) , (31) 
oe\ oe 2 oes 



df /e?Ve (see equation (18)) in terms of the order parameters reads: 



where e, is a canonical unit basis vector. The divergence of the higher-order stress term S = 

(32) 



V • S = V • {^ ^j e ^ ( e i ® ei - e 2 <8> e 2 ) = A e V 2 e 2 (ei ® ei - e 2 ® e 2 ) 



The chemical potential [i is of the form (see equation (|22[)): 

H = y c ~ AcV ' c ' ( 33 ) 

where df/dc is the 'local' part of the chemical potential. The time derivative of the entropy density, 
which will play a role in formulating the heat transport equation, reads: 

d 2 f ■ d 2 f . d 2 f . . d 2 f d 2 f „ 

s = —T — ■ k - -Ve — c - • Vr 

dT 2 dTde ' dTdVe' dTdc dTOVc ,„.x 

cy f d 2 f . d 2 f . d 2 f . W 

T rcf ar^ei 61 dTde 2 e2 dTdc^ 

where in the last line we have taken into account that temperature does not affect the phase 
boundary energy. 

5.3. Constitutive models for the dissipative terms 

Constitutive models for the mass flux, the heat flux and the viscous stress do not follow as a 
direct consequence of the chosen of the free-energy, but must be chosen such that the restrictions 
in Section [4] are satisfied. As is usual, the definition of a free-energy does not provide information 
on the kinetics (rate) of phase transformation processes. 

We opt for a viscous stress that depends on the rate of e 2 only, 

a v = ve 2 (ei (8) e\ — e.2 ® e 2 ) , (35) 



where v > 0, in which case equation (20 ) is satisfied. The mass flux j c is assumed to be proportional 
to the chemical potential, 

j c = -MVm, (36) 
where M is known as the 'mobility'. The expression 

M = M exp(-%), (37) 



T 
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where M$ > and Q > are constants, will be used in the simulations. The thermodynamic 
restriction in equation (23) is satisfied if M > 0. For the heat flux, the Fourier model is assumed: 

q = -kVT, (38) 



where k is the thermal conductivity. The thermodynamic restriction in equation (25) is satisfied 
if k > 0. 



6. Governing initial boundary value problems 



A complete set of governing equations for the considered model can now be defined. The model 
is based on momentum balance, mass balance for the solute phase and an energy balance for heat 
transport. The primal unknowns to be solved for are the displacement, the solute concentration 
and the temperature. 



The boundary dQ is assumed to be partitioned into subsets such that dT^ U dF g 



dT q U dT T 



dn and dT h n dT„ 



dT q n dT T 



The balance of momentum equation follows directly from 



equation Q. Neglecting differences in mass density between phases, the balance of momentum 
equation and the associated boundary and initial conditions read: 



pu — V • a 


= b 


in x [0, if] 


(39) 


an 


= h s 


on T h x [0, i f ] 


(40) 


u 


= 9u 


on F g x [0,i f ] 


(41) 


A e Ve2 • n 


= 


on (90 x [0, if] 


(42) 


u (x, 0) 


= u 


on O 


(43) 


ii (x, 0) 


= Vq 


on £1. 


(44) 



The boundary condition Ve2 • n = is motivated by the interpretation that the gradient terms 
provide a nonlocal effect, and in the absence of neighbouring material particles on the boundary 



no nonlocal work can be done on dQ.. This coined the 'insulation condition' by Polizzotto (2003). 



Inserting the expression for the mass flux in equation (36) into the differential expression of 
mass conservation §5§ leads to an equation governing the transport of the solute. For the form 
of the free-energy given in Section 5.1 and the definition of fi in equation (22) in terms of the 



free-energy, the mass diffusion equation and associated boundary and initial conditions read: 

r - V • \ [V ( ^ - AV 2 c J =0 on tl x [0, if] 



A c Vc • n = on dQ x [0, if] 
c (x, 0) = cq on f2. 



(45) 

(46) 

(47) 
(48) 



The above equation is known as the Cahn-Hilliard equation. A consequence of the manner in 
which the surface energy is included in the model is the presence of fourth-order derivatives. The 



boundary condition in (46) implies that there is no mass flux across dQ (see (22) and (36)), and 
(47) is interpreted in the same fashion as the condition on Ve2 • n in the momentum balance. A 



11 



more general form of the boundary conditions for the Cahn-Hilliard equation can be found in Wells 



et al. (2006). 



The final governing equation is the energy balance. Using the definition of the Helmholtz 
free-energy, u can be replaced by / + Ts + Ts in the energy balance (10), yielding 

/ + Ts + Ts = -V • q + a : Vv + (V • E) : £ + £ : Ve + (V • £)c + £ • Vc - /iV • j c - V// • j c - (49) 



By considering the expansion of / in (15) and the constitutive models for <r e , [i and s in equa- 
tions (19), (22) and (24), respectively, the energy balance reduces to 

Ts 



-V • q + a v : e - Vfi • j c . 



(50) 



Using the expression of s in equation (34), the heat equation, together with suitable boundary and 
initial conditions reads: 



Tc v 
T Te i 



T-T 



da . du 

— : e H — c 

dT dT 



+ V/i • j' c — crt, : e + V • q = on O x [0, if] 



(51) 



q-n = h T onT g x [0,t f ] (52) 
T = g T onr T x[0,t f ] (53) 
T(x,0)=T onfi. (54) 



The complete model involves the solution of all three coupled equations. 



7. Fully-discrete Galerkin formulation 

A Galerkin finite element formulation of the governing equations is now developed which will 
be used in computing numerical examples. A feature of the governing equations is the presence 
of fourth-order derivatives of the concentration and displacement fields. Ordinarily, this would 
require finding approximate solutions in a subspace of H 2 (fl), but such finite element element 
spaces are troublesome to construct. We will exploit two different strategies to avoid this difficulty. 
For the mass diffusion equation it is convenient to adopt a mixed formulation, which is essentially 
an operating-splitting approach and is not encumbered with the difficulties that plague mixed 
formulations that result in saddle-point problems. For the momentum balance, mixed schemes 
are difficult to analyse and in some cases stable schemes may not be known. To circumvent this 
difficulty, techniques for fourth-order elliptic equations that are inspired by discontinuous Galerkin 
methods will be used (see |Engel et al] ( |2002[ ); |Wells and Dung| fl2007| )). These methods permit 



the rigorous solution of fourth-order problems using a primal formulation and H 1 (fi)-conforming 
finite element spaces. 

To formulate a finite element problem, let T be a triangulation of f2 into finite element cells 
such that T = {K}. We will work with the usual Lagrange finite element basis 

V k = {v£ H 1 (O) , v £ P k (K)VK € T} , (55) 

where P k denotes the space of Lagrange polynomials of order A: on a finite element cell. The 
numerical formulation that we propose involves solving the following problem: given the data at 
time t n , find u h ,c h ,T h ,T h £ (V k2 ) d x V kl x V kl x V kl at time t n+ i such that 



L(w,q,r,z;u h ,c h ,T h ,T h ) =0 V (w, q, r, z) G (V k2 ) d x V kl x V kl x V kl , (56) 
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where k 2 > 1 and k\ > 0. The functional L is linear in each of (v,q,r,w) and therefore can be 
split additively into four contributions: 



L(w, q, r, z; u h , c h , r h , T h ) = L u (w; u h , c h , r h , T h ) + L c (q; u h , c h , r h , T h ) 

+ L T (r; u h , c h , T h ,T h ) + L T (z; u h , c h , T h ,T h ), (57) 

where L u , L c , L T and Lt will represent the contributions of linear momentum, mass diffusion, the 
chemical potential and heat transport, respectively. Each of these terms will be defined in this 
section. 

Implicit time integrators will be used for all equations. We believe this to be the only feasible 
approach owing to the parabolic nature of the mass and heat transport equations, and the presence 
of fourth-order terms in the momentum balance and mass diffusion equations. 



It is particularly convenient to express the problem using the format of equation (56) since 
we use high-level tools to generate computer code automatically for this problem, and these tools 
inherit this mathematical expressiveness. Furthermore, the tools perform automatic differentiation 
(both regular and directional derivatives) so there is no need to compute by hand a linearisation 
of the problem. The functional L will be the input to the code, and the directional derivative 
(the linearisation) is computed using automatic differentiation, yielding the Jacobian for use in a 
Newton solver. This will be expanded upon in the following section. 

The subscript '/i' will be used in this section to denote an approximate quantity, for example 
fh = f(uh), where Uh is the approximate displacement field. 

7.1. Balance of momentum 



Multiplying the balance of momentum equation ( 39 ) by a function w and applying integration 
by parts to various terms, 

pw ■ Uh dx + Vw : &h dx + / Vw : V • S^, dx + Vw : <r v ^ dx 
n Jn Jn Jn 



/ w bdx— / w-h s ds = 0, (58) 
Jn Jan 



where the decomposition of the stress into 'local' and 'gradient' contributions has been used, 
&he = — V • S/i, and the Neumann boundary condition an = h s has been inserted. Addressing 



now the term involving S/j and considering the precise form of 5]^ in (32), we have 



/ Vw : V--E h dx = - [ (w hl -w 2 ,2)K^ 2 e 2 ,hdx. (59) 
Jn Jn 

This term is problematic since after the application of integration by parts, second-order spatial 
derivatives of w and are present, which classically would require searching for approximate 
solutions using ^T 2 -conforming functions, thereby precluding the use of standard Lagrange finite 
element basis functions. To circumvent this difficulty while still using iT 1 -conforming element 
and without abandoning a primal approach, we introduce integrals over finite element cell facets 
to impose weak continuity of the normal derivative while preserving consistency and stability, as 



formulated for the biharmonic equation in Engel et al. (2002) and the Cahn-Hilliard equation in 
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Wells et al. (2006). This then permits the use of standard iJ 1 -conforming finite element basis 



functions. The modified variational form, using df/dei in place of <r, reads: 



pw ■ Uh,n+l-a m dx+ I (w 1} l + U7 2 , 2 ) 



df h ,n+l— at 



+ / (m,2 + w 2 ,i) 



df h , 



n 

n+l— a f 



dx + 



, x 9fh,n+l— at 
[Wn - W 22 ) q dx 

oe 2 



de 3 



dx + y I ■ 

IK 



V (wi,i - w 2)2 ) ■ A e Ve 2 ,h,n+i-a f dx 



Jjw 1A - w 2j2 j • (K^e 2) h,n+i- af ) ds - j[(\ e V (wi 7 i - w 2}2 )) ■ le 2 ,h,n+i-a f j ds 
+ j_ ^ [^1,1 - u> 2 ,2l • W,h,n+l-af\ dT + J (wx t i - w 2 , 2 )ve 2Xn j rl _ af dx 



w ■ b n+ i_ a dx 



ds, (60) 

Jdn 



where T denotes the set of all interior facets, [a] = a+n + + a__n_ and (b) = (5+ + b ) /2, '+' and 
' — ' denote opposite sides of a cell facet, rj is a dimensionless penalty parameter and /i^ is a measure 
of the local element size. The penalty parameter is required for stability of the formulation and 
is of order one. Time derivatives are dealt with using a generalised-a scheme (see, for example, 



Chung and Hulbert (1993)). The acceleration, velocity and displacement at the end of a time step 

1 



are related via 



u n+ \ = u n + Atu n + 



P )u n + /3ii„+i 



u n+1 = u n + At ((1 - 7) u n + ju n+1 ) , 



(61) 
(62) 



where At = £ n+ i —t n , (3 and 7 are parameters, and mid-point values of primal fields are computed 
according to 

Vn+i-a = (1 - a) y n+ i + ay n . (63) 

where a is a parameter. The parameters a m , a.f, 7 and /3 determine properties of the time stepping 
scheme and will be reported with other problem data for the numerical examples in the following 
section. Nonlinear functions are evaluated at a mid-point, for example 



n+a ) ■ 



(64) 



The discontinuous Galerkin-type approach is consistent for k > 1, and analysis of the bihar- 
monic equation has shown that the method is stable for sufficiently large a (typically a is of 
order 1). It has been proven that the method converges in the L 2 for the biharmonic equation at 
a rate of k + 1 for k > 2 ( |Engel et alj |2002[ ) and at a rate of k for k = 2 ( |Wells and Dung[ |2007j ). 

7.2. Mixed form of the Cahn-Hilliard equation 

For the Cahn-Hilliard equation we adopt an operator splitting approach to deal with the fourth- 
order spatial derivative and split equation (45) into two second-order equations: 

c - V ■ MVr = 0, 



% + A c V 2 c = 0, 
oc 



(65) 
(66) 
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where c and r are the independent unknowns that will be solved for. If the above problem is solved 



exactly, then r and [i will coincide (see equation (22)). It will be important to distinguish between 
r and /U when finding approximate solutions. Note that through df /dc there is a dependency on 
the strain and the temperature. 

Casting the above set of equations into a weak form, using the Crank-Nicolson method for the 



time derivative, and applying the boundary conditions from equations (47) and (46 ), the functionals 
L c and L T read: 



Q Ch,n+ \, Ck,n dx + / Vg • MVr M+1 /2 dx, 
q At Jn 



(67) 



df, 



Lfj,= I rT hjTl+1 dx - [ r a/ ^ n+1 dx - I Vr • A c Vc hjTl+ i dx 
Jn Jn ° c Jn 

This operator splitting approach to the Cahn-Hilliard equation was presented and analysed by 



Elliott et al. (1989), and for the form of the free-energy and mobility considered by Elliott et al. 



(1989) it was shown to be stable. A discontinuous Galerkin type approach, similar to that used for 



the momentum equation, can also be used to solve the Cahn-Hilliard equation in its primal form 



(Wells et al., 2006). 



7.3. Heat transport equation 



The heat equation (51) can be cast into a weak form via the usual process. Using the Crank- 



Nicolson method to deal with time derivatives, the functional Lt reads: 

Th n.4-1 /•}. Th n.4-1 — Th n 



h,n+l/2 J-h,n+l 
< > T re f At 



dx 



zT, 



h,n+l/2 



9cr e ,h,n+l/2 _ £h,n+l - £h,n . Ch, n +l ~ Ch 



dT At 



+ In Z ' Vr?i '" +1/2 ' ■ 7 '™+ 1 /2 



+ 



r D,n+l/2 



dT At 



dx 



At 



dx 



/ Vz-q h >n+ i/ 2 dx- / zh q , 2 ds. (69) 



Taking into account details of the constitutive models, 

3 



Th,n+l/2 Th, n +1 — Th,n 



n 



! ref 



At 



dx -^T/,, n+ i/2 



i=l 



9 f &i,h,n+l &i,h,n 

deST At 



dx 



z Ti 



9 2 f Ci hn+l ~ Cih; 



o v h ' n+1/2 dc8T 

zh>e 2,h,n+l/2 



At 

e-2,h,n+l — Z2,h,n 



+ Vr M+1/2 • M h VT hin+1/2 dx 



At 



dx+ I Vz- k T VT h n+1 , 2 dx 
n 



Z K,n + l/2 dS - ( 70 ) 



8. Numerical examples 

A number of numerical examples are now presented to demonstrate that the model can quali- 
tatively capture classical processes observed in the heat treatment of steels. It is not the intention 
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Figure 3: Mobility as a function of temperature. 

to consider realistic material parameters at this stage. The determination of physical parameters, 
and computing with these, is a substantial undertaking and is the subject of ongoing work. 

The model includes a number of coupled processes, which permits considerable freedom in 
how transformations are induced, and a variety of parameters play a role in the development of 
microstructure. The presented examples involve varying the heat flux across the boundary of the 
domain for a fixed set of model parameters. For the parameters appearing in the free-energy, the 
following values are used: A A = 7.3 x 1(T 3 ; A 2 = 6.6 x 1(T 3 ; T P = 1; A c = 1 x 1(T 7 ; B 62 = 0.3; 
B 42 = 3.1 x 1(T 3 ; B 22 = 2.5 x 1(T 3 ; B 3 = 5 x 1(T 2 ; A e = 1 x 1(T 8 ; T M = 0.495; a = 1 x 1(T 2 ; 
x\ c = 5 x 1CT 3 ; x\2 = 1 x 10 _1 ; ct = 3.5 x 1CT 8 ; and x 2c = 8 x 10 -2 . For model parameters 
that do not appear in the free-energy, the following values are used: p = 5 x 10~ 7 ; !/ = 5x 10~ 8 ; 
kr = 2 x 10~ 2 ; Mq = 1 x 10 4 ; and Q = 5. The adopted time stepping parameters are p^ = 0.7; 
a m = (2poo - l)/(poo + 1); a f = Poo/(Poo + 1); /? = (1 - a m + «/) 2 /4; and 7 = 1/2 - a m + a f . 
The time stepping parameters are chosen such that the momentum balance scheme is second- 
order accurate and strongly stable (at least for linear problems). The mobility M as a function of 
temperature for the adopted parameters is plotted in Figure |3j There is considerable scope for a 
extensive computational studies into the impact of various parameters on different microstructural 
processes. 

All simulations are performed on a unit domain £1 = (0, 1) x (0, 1) using triangular finite element 
cells in a regular diagonal pattern and with 128 vertices in each direction. For the displacement 
field, T g = (x, 0) U (0, y) £ dtt with it = on T g and h s = on Y h = dVL\T g . With this boundary 
condition, any effect of a displacement constraint on microstructure evolution will be evident. The 
heat flux across dVt, which we denoted by h q , is set proportional to the difference between the 
temperature T on the boundary and a prescribed 'external' temperature T ex t: 

hq = —S (T — Text) , (71) 

where 5 > is a parameter. The larger the value of 6, the more rapid the cooling. Other 
boundary conditions were defined in Section[6} The impact of the cooling rate on the microstructure 
development will be examined by varying 5 and T ex t • The initial conditions for the displacement and 
concentration are fields that are perturbed randomly about zero. In the case of the concentration 
field, Co £ [— 10 -4 , 10 -4 ] is a random number with a uniform distribution, and for the initial 

displacement, uq G [— 10~ 4 /i^, 10 -4 /i^-] 2 is a random vector with a uniform distribution and hj( 
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Figure 4: Contours of the displacive order parameter for the rapid cooling case to T GXt =0.1. The approximate 
range of temperatures in the domain at the snapshots are: at t = 0.01, T ~ 0.27 — 0.61; at t = 0.015, T ~ 0.17 — 0.31; 
at t = 0.025, T w 0.1; and at t = 0.3, T « 0.1. 



is a measure of the finite element cell size. The dependency on the cell size is introduced so 
that the magnitude of the initial strains is roughly equal on different meshes. In all simulations 
T rcf = T = 1.2. 

Quadratic Lagrange functions are used for the displacement fielcQand linear Lagrange functions 
are used for all other fields (&2 = 2 and k\ = 1 in equation (56)). A fully-coupled solution strategy 
is used and a Newton-Krylov method is employed to solve the nonlinear equations in each time 
step. For the stabilising penalty term, r\ = 8. The problem-specific parts of the computer code 
used to perform the simulations have been generated automatically from a high-level description 
that resembles closely the notation used in this work by using a number of tools from the FEniCS 



Project (Logg and Wells, 2010 01gaard and Wells, 2010 01gaard et al. , 2008). The complete 



computer code used to perform all simulations reported in this work consists of one file only and is 



freely available under a GNU public license for both scrutiny and use (Wells and Maraldi 2011) 



8.1. Rapid boundary cooling to T ext = 0.1 

We first consider rapid cooling from an initial uniform temperature of T = 1.2, at which 
temperature austenite is the stable phase, down to T ext = 0.1, which is well-below the martensitic 
transition temperature of Tm = 0.495. At T = 0.1, the mobility of the diffusive phase is negligible, 
hence for this rapid cooling case the formation of a diffusive phase is not anticipated; the resulting 
structure will therefore be martensite. To simulate rapid cooling, the boundary heat flux parameter 
5 = 5 x 10 -2 is used. For this case, simulations were performed with an initial time step of 
At = 1 x 10~ 4 , which was increased to At = 2 x 10 3 at t = 0.266. The computed order parameter 
e2 at various times during the cooling process is shown in Figure [4j The diffusive variable c (not 
shown) is very close to zero everywhere in the domain. The rapid formation of fine martensite 
twins can be observed in Figure |4j Note that the martensite twins are slightly less developed at 
the left-hand and lower edges of the domain, due to the constraint u = 0, which inhibits lattice 
distortion on these constrained boundaries. This is also what induces the bottom left-hand corner 
to top right-hand corner orientation of the twins, rather than at 90° to this orientation. 



2 It is known that the formulation we have adopted to deal with the fourth-order derivatives of the displacement 
field leads to order two convergence in the L 2 -norm for the biharmonic equation when fs2 = 2 (Wells and Dung 



20071, which may appear to be a sub-optimal choice. However, we are interested primarily in e-2, and the adopted 



method converges with order two in the i/ 1 -norm when k-z = 2. 
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i = 0.0175 t = 0.02 t = 0.05 t = 0.5 



Figure 5: Contours of the displacive order parameter e2 for the rapid cooling case (5 = 5 X 1CT 2 ) to T cxt = 0.3. The 
approximate range of temperatures in the domain at the snapshots are: at t = 0.0175, T ~ 0.35 — 0.45; at t — 0.02, 
T » 0.32 - 0.41; at t = 0.05, T w 0.3; and at £ = 0.5, T » 0.3. 

S.S. Boundary cooling to T ex t = 0.3 ai different rates 

We now consider a less deep quench to T ex t = 0.3, which is still below the martensitic transition 
temperature, with the same boundary heat flux parameter of 5 = 5 x 10~ 2 . At T = 0.3, the mobility 
of the diffusive phase is still negligible. For this case, simulations were performed with a time step 
of At = 2.5 x 10 -4 . The order parameter e2 at various times during the cooling process is shown in 
Figure[5j For this case, the field c (not shown) is close to zero throughout the simulation, despite the 
diffusive phase being energetically favourable over the displacive phase, because diffusive processes 
are unable to develop due to the brief length of time over which the temperature, and hence the 
mobility, is sufficiently high. As with the T ex t = 0.1 case, martensite forms quickly, although it 
is now slightly coarser in the early stages. On a slightly longer time scale, there is some limited 
twin boundary motion as part of a minimal coarsening process. Such coarsening is not observed 
at lower temperatures, and can be attributed to the different balance between bulk and surface 
energy (in the model, the surface energy is unaffected by temperature while the bulk energy is). 
Noteworthy is that the microstructure that forms in this case is more affected by the displacement 
constraint on the two sides of the domain. This can be explained by the smaller thermodynamic 
force driving the e^ patterning compared to the T ext = 0.1 case. 

For this cooling case, the temperature contours are shown at two snapshots in Figure |6j In the 
absence of thermal coupling effects, the temperature contours would be symmetric about the x- 
and y-axes. Figure [6] shows a subtle breaking of symmetry. In particular, at t = 0.02 an alignment 
of contours with the direction of the martensite twins can be detected. The presented simulations 
are driven by large changes in the boundary heat flux, which makes it difficult to detect local 
temperature changes due to phase changes since these are relatively small, but also with possibly 
high gradients that are smoothed rapidly by conduction. The heat capacity plays a large role in 
temperature changes since it determines the change in temperature associated with a given energy 
release during a transformation. 

An intermediate cooling rate case with 5 = 1.5 X 10 -3 and a time step of At = 2 x 10~ 3 is 
now considered. For this case, the order parameters ei and c are shown in Figure [7] at a number if 
time steps. At this cooling rate, the impact of non-isothermal effects are clear. There is sufficient 
time for a diffusive phase to develop in three of the corners. The formation of a diffusive phase in 
the corners is accelerated by the combination thermoelastic effects and the boundary constraints 
and the temperature changes. Upon cooling, regions of volumetric straining with opposite signs 
develop in three of the corners which, due to the x\ c ^ 0, encourages diffusion of c. However, 
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t = 0.0175 t = 0.02 

Figure 6: Temperature contours for the rapid cooling case (5 = 5 x 1(T 2 ) to T cxt = 0.3. Thermal coupling effects 
lead to the subtle breaking of symmetry in the temperature contours. Symmetry is restored at later times via heat 
conduction. 



before a diffusive phase can form in the entire domain, the temperature drops to a level at which 
the mobility is too low for diffusion to continue, and martensite twins form in the regions in which 
the diffusive phase has not yet developed. For t > 1, the phase configuration is essentially stable 
and no further changes could be detected. Figure [8] shows the temperature contours for this case 
at four time instants. Initially a bias can be detected with the warmest region running between the 
two regions with the most developed diffusive phases (at t = 0.25). As time progresses, a subtle 
alignment of the temperature contours can be detected. At long times, the temperature field is 
smoothed by conduction, but at t = 0.2 very subtle temperature variations can be detected, with 
the temperature contours aligned with the martensite twins. This can be attributed to very slow 
twin coarsening effects. 

Even slower cooling rate cases are now considered, from which the impact of the cooling rate 
on the structure of the diffusive phase can be observed. For 5 = 1 x 10~ 3 and At = 2 x 10 -3 , 
the evolution of the diffusive phase is shown in Figure [9j At this cooling rate, a reasonably fine 
microstructure develops, with considerable domain boundary coherence. That is, the lamella near 
the boundaries are aligned with the boundary. Between t = 0.3 and t = 0.4 the temperature in the 
domain increases due to the energy released during the transformation, and the time is insufficient 
for conduction to transport the heat across the boundary of the domain. The microstructural 



features change if the cooling rate is slowed further, as shown in Figure 10, where S = 1 x 10~ 4 
and At = 1 x 10 -2 . The microstructure is clearly coarser, and the boundary coherence is reduced. 
Beyond t = 8, further microstructural changes could not be perceived. 

8.3. Rapid boundary cooling to T ext = 0.4 

The domain is now cooled to T = 0.4 with 5 = 5 x 10 -2 . The time step was changed during the 
simulation, and rang ed between At = 2.5 x 10" 4 and At = 5 x 10~ 2 . At T = 0.4, diffusive processes 
are very slow, but not yet negligible. For this case, the computed contours of and c are presented 



in Figure 11 It is clear that martensite twins form initially. The thermodynamic driving force for 
the development of a diffusive phase is greater than for the displacive phase, although it is inhibited 
by the low mobility, but eventually diffusive process can be detected and the diffusive phase slowly 
replaces the displacive phase. At this external temperature, the structure of the diffusive phase is 
quite fine. 
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t = 0.4 



t = 0.5 



t = 2 



Figure 7: The displacive order parameter ei and the diffusive order parameter c for the moderate cooling rate case 
to Text = 0.3. The approximate temperature in the domain at the snapshots is: at t = 0.25, T ~ 0.62; at t = 0.4, 
T w 0.45: at t = 0.5, T w 0.38: and at t = 2, T w 0.3. 





t = 0.4 




Figure 8: Temperature contours for the moderate cooling rate case to T ext = 0.3. Variations of the temperature con- 
tours away from a symmetric pattern are attributable to thermal coupling effects. At t = 2, very subtle temperature 
contours aligned with the martensite twins can be detected. 
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Figure 9: Contours of the diffusive order parameter c for the slow cooling case (<5 = 1 X 1CT 3 ) to T cxt = 0.3. The 
approximate temperature in the domain at the snapshots is: at t — 0.3, T w 0.72; at t = 0.4, T « 0.76; at t = 0.6, 
T w 0.66; and at t = 4, T w 0.3. 




0.6 



-0.6 



t = 1.75 



t = 2.5 



t = 4 



t = 8 



Figure 10: Contours of the diffusive order parameter c for the very slow cooling case (5 = 1 X 10 4 ) to T cx t = 0.3. 
The approximate temperature in the domain at the snapshots is: at t = 1.75, T w 0.91; at t = 2.5, T w 0.89; at 
t = 4, T « 0.8; and at * = 8, T w 0.6. 
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Figure 11: The displacive order parameter e2 and the diffusive order parameter c for the rapid cooling case to 
Text = 0.4. At all snapshots, the temperature in the domain is approximately 0.4. 
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Figure 12: The displacive order parameter e2 and the diffusive order parameter c for the rapid cooling case (<5 
5 x 10~ 2 ) to T cx t = 0.45. At all snapshots, the temperature in the domain is approximately 0.45. 



An extremely challenging aspect of computations at this temperature is reconciling the dramatic 
time scale differences for the different processes, especially during the latter stages. The diffusive 
phase develops very slowly, but the evolution of the diffusive phase can induce sporadic changes in 
the martensitic phase. Using a time step that is suitable for diffusive processes can lead to a loss 
of solution stability when sporadic displacive changes take place. There is considerable scope for 
developing efficient, robust and accurate methods for spanning the large difference in time scales. 



8-4- Rapid boundary cooling to T ext = 0.45 
The final example is rapid cooling with 5 



5 x 10 2 to T cxt = 0.45. The time step was changed 



during the simulation, and ranged between At = 5 x 10 4 and At = 2 x 10 1 . Snapshots of the 



order parameters e 2 and c are presented in Figure 12 As expected, martensite develops quickly. 
Interestingly, the domain then undergoes a rearrangement between t = 0.1 and t = 0.4, in which 
e 2 goes from an alternating pattern to a value close to zero in the bottom left-hand region of the 
domain. It is likely that this is due to the boundary constraint resisting deformation, and the 
driving force for the square-to-rectangular transition being too small to overcome this constraint. 
In the top right-hand triangle, e 2 is approximately constant at t = 0.4. At T = 0.45, the driving 
force behind the development of the alternating e 2 pattern is relatively weak, whereas the energy 
associated with phase boundaries does not have a direct temperature dependency. The system 
undergoes another change between t = 0.4 and t = 2, with the sign of e 2 in the upper right- 
hand region changing. On a longer time scale, diffusive processes develop, first in the region where 
e 2 ~ and then progressing in into the region where e 2 is not equal to zero, eventually replacing the 
martensite phase. Noteworthy is the lack of structure in the diffusive phase in the bottom left-hand 
triangular region, and the more coherent lamella structure in the upper right-hand region. 

As for the previous example, computing solutions to this problem robustly and on a reasonable 
computational time scale is extremely challenging. 



9. Conclusions 

A phase field model that can simulate both diffusive and displacive phase transitions has been 
presented. The model is developed in a formal thermodynamic setting, with free-energy expres- 
sions associated with distinct microstructural processes, the postulation of an energy balance and 

22 



demonstration that the classical entropy inequality is satisfied. Various couplings terms introduce 
dependencies between different processes. The model includes both bulk energy and phase bound- 
ary energy, which impacts on the fineness of the computed microstructures. Surface energies have 
been introduced via higher-order spatial gradients of the relevant order parameters, and this de- 
mands special care in the formulation of the thermodynamic balance laws. To satisfy the classical 
point-wise form of the entropy inequality, non-standard stress-like terms have been introduced to 
the energy balance equation. These non-local terms are related to phase boundary surface energy, 
and represent the work done by particles that neighbour an arbitrary subdomain. 

The differential equations that follow from the fundamental balance laws involve the displace- 
ment field, the solute concentration and the temperature field, with significant couplings between 
all equations. The presence of surface energies in the model leads to a momentum balance and 
mass balance (diffusion) equations that involve fourth-order spatial derivatives. The presence of 
the higher-order derivatives is problematic when developing a Galerkin finite element method for 
solving the problem. To counter this, a sophisticated Galerkin method has been formulated for the 
problem that imposes the classically required solution regularity in a weak sense. The computer 
model has been generated for a large part automatically from a high-level specification, and is 
published as supporting material. 

It has been demonstrated through numerical simulations that the proposed model can capture 
qualitatively a variety of observed phenomena, including the formation of martensite twins, the 
development of pearlitic structures and pearlite replacing martensite when a specimen is held at 
a sufficiently high temperature for an extended period. Phase transitions have been triggered via 
control of heat flux across the boundary, but transformations can also be induced by mechanical 
loading. The application of the model with physically determined parameters requires further 
investigation. Such an investigation will also necessitate further development of numerical solution 
strategies to enable simulations on large domains and appropriate time scales to be performed both 
accurately and within a tolerable simulation time. 
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